##Replication Code- Michael O'Neill and Amy Uden; March 2016

##Replicating "Polarization and the Decline of the American Floating Voter"
##by Corwin Smidt

#####################################################
#####################################################

#Replication Code for Figures 1-3:

##Replication Code for Figures 1-3
##Note: This code closely follows Smidt's code in order to ensure
##visual comparability of replication figures. I walked through each
##step of the code to see whether it made sense (in .do files and graphing code for R). 
##I also checked to ensure figures replicated.

setwd(dirname(file.choose()))

require(foreign)

data1 <- read.csv("modifiedcdf.csv")

head(data1)
nrow(data1)

summary(data1)

##Replication Code for figures:
##accepting Smidt's compiled dataset, the following figures replicate:

##Figure 1:
##First, I need to run the create.do file. I have reviewed Smidt's code 
##and rerun the figure1_create.do file as my own .do file ("replication_figure1_create.do"). 
##The output file, which includes data needed to create figure 1, is saved as figure1.dta.
##I used stata version 12 because I have found it more readily compatible with R. 

##Now, developing my code from Smidt's code:

rm(list=ls())

setwd(dirname(file.choose()))

library(foreign)
library(ggplot2)

statadata <- read.dta("figure1.dta")

catdata <- data.frame(Year=statadata$year,Percentage=statadata$percentage,Behavior=statadata$category)

plot <- qplot(Year,Percentage,data=catdata,geom="line",linetype=Behavior) + annotate("text",label="Standpatters", x=1970,y=57,size=4)  + annotate("text",label="Floating Voters", x=1961,y=6,size=4)  + annotate("text",label="Repeat Non-Voters", x=1989,y=27,size=4) + annotate("text",label="Surge & Decliners", x=1963,y=24,size=4)

pdf("figure1.pdf", width=8,height=6)
plot + theme_bw() + scale_y_continuous(limits=c(5,60),breaks=c(10,20,30,40,50,60)) + labs(linetype="Voting Behavior")  + theme(legend.position="none")
dev.off()

####################################

####Figure 2:

statadata2 <- read.dta("figure2a.dta")

sophplotdata <- data.frame(Year=rep(statadata2$year,2),Group=c(rep("Low", 15),rep("High", 15)), factor = rep("By Political Awareness",30), Percent = c(statadata2$low,statadata2$high))

pidplotdata <- data.frame(Year=rep(statadata2$year,3),Group=c(rep("Independent",15),rep("Lean/Weak",15),rep("Strong",15)), factor = rep("By Strength of Partisanship",45), Percent = c(statadata2$ind, statadata2$weak,statadata2$strong))

sophplotdata$Percent <- sophplotdata$Percent*100
pidplotdata$Percent <- pidplotdata$Percent*100

combinedata <- rbind(sophplotdata,pidplotdata)


lp1 <- c("High"=1,"Low"=2)
lp2 <- c("Strong"=1,"Lean/Weak"=2,"Independent"=3)


sophplot <- qplot(Year,Percent,data=sophplotdata,geom="line",linetype=Group,xlab="Year", ylab="Percentage")  + scale_x_continuous(limits=c(1960,2012),breaks=c(1960,1968,1976,1984,1992,2000,2008)) + scale_y_continuous(limits=c(20,95),breaks=c(20,40,60,80)) + scale_linetype_manual(name="By Political Awareness:",values=lp1) + theme_bw() + theme(legend.position=c(.3,.85),legend.background = element_rect(colour = "grey"))  + theme(axis.title.y=element_blank(),axis.text.y=element_blank(),axis.ticks.y=element_blank())
pidplot <- qplot(Year,Percent,data=pidplotdata,geom="line",linetype=Group,xlab="Year", ylab="Percentage") + scale_x_continuous(limits=c(1960,2012),breaks=c(1960,1968,1976,1984,1992,2000,2008)) + scale_y_continuous(limits=c(20,95),breaks=c(20,40,60,80)) + scale_linetype_manual(name="By Strength of Partisanship:",values=lp2) + theme_bw() + theme(legend.position=c(.33,.89),legend.background = element_rect(colour = "grey")) 
#bothplot <- qplot(Year,Percent,data=combinedata,geom="line",linetype=Group,colour=factor,xlab="Year", ylab="Percentage") + scale_x_continuous(breaks=c(1952,1960,1968,1976,1984,1992,2000,2008)) + facet_grid(.~factor) + theme_bw() + scale_linetype_manual(values=lp1)  + scale_linetype_manual(values=lp2)


library(grid)
pdf("figure2.pdf", width=8,height=6)
grid.newpage()
pushViewport(viewport(layout=grid.layout(1,100)))
vplayout <- function(x,y)
  viewport(layout.pos.row=x,layout.pos.col=y)
print(pidplot,vp=vplayout(1,1:51))
print(sophplot,vp=vplayout(1,52:100))
dev.off()

####################################

####Figure 3:

issuereg <- read.dta("figure3a.dta")

sophplotdata <- data.frame(Year=rep(issuereg$year[4:15],2),Group=c(rep("Fairly Low", 12),rep("Fairly High", 12)), factor = rep("By Political Awareness",24), Awareness = c(issuereg$low[4:15],issuereg$high[4:15]),lb = c(issuereg$lowlb[4:15],issuereg$highlb[4:15]), ub=c(issuereg$lowub[4:15],issuereg$highub[4:15]))
pidplotdata <- data.frame(Year=rep(issuereg$year,2),Group=c(rep("Independent",15),rep("Strong",15)), factor = rep("By Strength of Partisanship",30), Awareness = c(issuereg$ind,issuereg$spid),lb = c(issuereg$indlb,issuereg$spidlb), ub=c(issuereg$indub,issuereg$spidub))


lp1 <- c("Fairly High"=1,"Fairly Low"=2)
lp2 <- c("Strong"=1,"Independent"=2)


sophplot <- qplot(Year,Awareness,data=sophplotdata,geom="line",linetype=Group,xlab="Year", ylab="Issue Awareness")  + scale_x_continuous(limits=c(1956,2012),breaks=c(1956,1964,1972,1980,1988,1996,2004,2012)) + scale_y_continuous(limits=c(20,90),breaks=c(20,30,40,50,60,70,80,90)) + scale_linetype_manual(name="By Political Awareness:",values=lp1) + theme_bw() + theme(legend.position=c(.3,.85),legend.background = element_rect(colour = "grey"))
#sophplot <- ggplot(sophplotdata,aes(Year,Awareness,ymin=lb,ymax=ub)) + geom_line(aes(linetype=Group))  + scale_x_continuous(limits=c(1956,2012),breaks=c(1956,1964,1972,1980,1988,1996,2004,2012)) + scale_y_continuous(limits=c(20,90),breaks=c(20,30,40,50,60,70,80,90)) + scale_linetype_manual(name="By Political Awareness:",values=lp1) + theme_bw() + theme(legend.position=c(.3,.85),legend.background = element_rect(colour = "grey"))
limits <- aes(ymax=ub,ymin=lb)
sophplot <- sophplot + geom_ribbon(limits,colour="lightgrey",alpha=.1,show_guide=FALSE)  + theme(axis.title.y=element_blank(),axis.text.y=element_blank(),axis.ticks.y=element_blank())

pidplot <- qplot(Year,Awareness,data=pidplotdata,geom="line",linetype=Group,xlab="Year", ylab="Awareness") + scale_x_continuous(limits=c(1956,2012),breaks=c(1956,1964,1972,1980,1988,1996,2004,2012))  + scale_y_continuous(limits=c(20,90),breaks=c(20,30,40,50,60,70,80,90)) + scale_linetype_manual(name="By Strength of Partisanship:",values=lp2) + theme_bw() + theme(legend.position=c(.35,.85),legend.background = element_rect(colour = "grey"))
limits <- aes(ymax=ub,ymin=lb)
pidplot <- pidplot + geom_ribbon(limits,colour="lightgrey",alpha=.1,show_guide=FALSE)

library(grid)
pdf("figure3.pdf", width=8,height=6)
grid.newpage()
pushViewport(viewport(layout=grid.layout(1,50)))
vplayout <- function(x,y)
  viewport(layout.pos.row=x,layout.pos.col=y)
print(pidplot,vp=vplayout(1,1:26))
print(sophplot,vp=vplayout(1,27:50))
dev.off()

######################################################
######################################################

##Replication Code for Table 1.

load("modifiedcdf.RData")

##or:

x <- read.dta("modifiedcdf.dta")

summary(x)
head(x)

##loading various packages and dependencies for bivariate probit
source("https://bioconductor.org/biocLite.R")
biocLite()
biocLite("graph")
biocLite("Rgraphviz")
install.packages("Zelig")
install.packages("ZeligChoice")
library("Zelig")
library("ZeligChoice")
library("foreign")


#################################################
#prepare data for fitting bivariate probit model

x<-cbind.data.frame(x["ambivalent"],x["undecided"],x["education"],x["age"],x["female"],x["sumissue"],x["black"],x["totalmentions"],x["pidstrength2"],x["changerdpi"],x["reelectcamp"],x["year"],x["weight"])

#check if recurisive and drop missing data
is.recursive(x)          
y<-as.data.frame(x[complete.cases(x),]) 

##checking if this is correct form###
summary(y)               
head(y)     

##Fit model with no weighting to replicate
out<- zelig(cbind(y$ambivalent, y$undecided) ~ y$education + y$age + y$female + y$sumissue + y$black + y$totalmentions + y$pidstrength2 + y$changerdpi +y$reelectcamp,
            model = "bprobit", data = y)

#initialize variables to bootstrap coefficients and standard errors
bootdat <- NULL
onecluster <- NULL
indexdat <- NULL
reps <- 2000
coeffs  <- matrix(NA, nrow = reps, ncol = 21)

set.seed(1004)

clusterboot <- for(i in 1:reps){
  
  bootdat <- NULL  
  model.sim <- NULL
  
  ##pull clusters by years with replace
  clusters <- names(table(y$year))
  draws <- sample( 1:length(clusters), length(clusters), replace=TRUE)  
  inds <- clusters[draws]
  print(i)
  
  
  for(j in 1:length(inds)){
    onecluster <- y[as.character(y$year) == inds[j], ]
    
    bootdat <- rbind(bootdat, onecluster)
  }
  
  
  model.sim <- zelig(cbind(bootdat$ambivalent, bootdat$undecided) ~ bootdat$education + bootdat$age + bootdat$female + bootdat$sumissue + bootdat$black + bootdat$totalmentions + bootdat$pidstrength2 + bootdat$changerdpi +bootdat$reelectcamp,
                     model = "bprobit", data = bootdat)
  
  coeffs.long <- as.data.frame(model.sim$getcoef())[ ,1]
  
  coeffs[i, ] <- t(as.matrix(coeffs.long[1:21]))
  
}

#I should save this considering how long that took...
save(coeffs, file = "BC_Coeffs_Table1_32116.Rdata")

##get mean of coeffs
e.coeff<-NULL
for(i in 1:ncol(coeffs)){ 
  v <-coeffs[, i]  
  e.coeff[i] <- mean(v)
}

e.coeff

##get standard errors from bootstrap
se <- NULL

for(i in 1:ncol(coeffs)){ 
  vector <-coeffs[, i]  
  se[i] <- sd(vector)
}

se


#p-values
t<-NULL
pv<-NULL

for(i in 1:21){
  mod <- e.coeff[i]
  mod.se <-  se[i]
  t <- (mod-0)/mod.se
  pv[i] <- 2*(1 - pnorm(abs(t), 0, 1))
}

pv

#bias
#Stata nebulously outputs a bias correction, possibly for bootstrapping.
#Note that this correction is not to be subtracted as the author does in the Stata help files.
#As this is a replication, the correction is here.

bias.correction<-NULL
bias.correction<-as.data.frame(out$getcoef())[,1]-t(as.data.frame(e.coeff))
bc.coeff<-e.coeff+bias.correction
t.bias.correction <- t(bias.correction)

c.coeffs<-coeffs

for(i in 1:21){
  c.coeffs[i, ] <-  coeffs[i, ] + bias.correction
}
head(c.coeffs)
#previous section with corrections;

bc.coeff<-NULL
for(i in 1:ncol(c.coeffs)){ 
  v <-c.coeffs[, i]  
  bc.coeff[i] <- mean(v)
}

bc.coeff



##get standard errors from bootstrap
bc.se<-NULL
for(i in 1:ncol(c.coeffs)){ 
  vector <-c.coeffs[, i]  
  bc.se[i] <- sd(vector)
}

bc.se


#p-values
bc.pv<-NULL
for(i in 1:21){
  bc.mod <- c.coeffs[i]
  bc.mod.se <-  bc.se[i]
  bc.t <- (bc.mod-0)/bc.mod.se
  bc.pv[i] <- 2*(1 - pnorm(abs(bc.t), 0, 1))
}

bc.pv

#Calculate aROC for each group


###############################################
##cross-era effects#attempt at

#These means are pulled from the the source Stata code and cannot be replicated through any permutation of mean or weighted mean, >= or >, or grouping by ambivalent v. undecided.
#Using these values against the 2000 and above values they should be of pre 1980 existant in Stata code with no source
#Given the nebulous origin of the cut-off values and ad hoc rescaling to normal distribution of a non-normal DGP this is a prime canidate for extension.

aeducation<-0
aage<-0
afem<-0
asumissue<-0
ablack<-0
atotalmentions<-0
apidstrength2<-0
achangerdpi<-0
areelectcamp<-0

for(i in 1:2000){
  
  
  
  yed80<-y
  yed80$education[y$year >1980]<-1.08453
  mi<-min(c.coeffs[i,1]+c.coeffs[i,4]*yed80$education + c.coeffs[i,6]*y$age + c.coeffs[i,8]*y$female + c.coeffs[i,10]*y$sumissue + c.coeffs[i,12]*y$black + c.coeffs[i,14]*y$totalmentions + c.coeffs[i,16]*y$pidstrength2  + c.coeffs[i,18]*y$changerdpi + c.coeffs[i,20]*y$reelectcamp)
  yed80$education[y$year >1980]<-1.762244
  ma<-max(c.coeffs[i,1]+c.coeffs[i,4]*yed80$education + c.coeffs[i,6]*y$age + c.coeffs[i,8]*y$female + c.coeffs[i,10]*y$sumissue + c.coeffs[i,12]*y$black + c.coeffs[i,14]*y$totalmentions + c.coeffs[i,16]*y$pidstrength2  + c.coeffs[i,18]*y$changerdpi + c.coeffs[i,20]*y$reelectcamp)
  aeducation[y$year[i]>=2000]<-aeducation+ma-mi
  
  
  yed80<-y
  yed80$age[y$year >1980]<-45.37945
  mi<-min(c.coeffs[i,1]+c.coeffs[i,4]*y$education + c.coeffs[i,6]*yed80$age + c.coeffs[i,8]*y$female + c.coeffs[i,10]*y$sumissue + c.coeffs[i,12]*y$black + c.coeffs[i,14]*y$totalmentions + c.coeffs[i,16]*y$pidstrength2  + c.coeffs[i,18]*y$changerdpi + c.coeffs[i,20]*y$reelectcamp)
  yed80$age[y$year >1980]<-46.31365
  ma<-max(c.coeffs[i,1]+c.coeffs[i,4]*y$education + c.coeffs[i,6]*yed80$age + c.coeffs[i,8]*y$female + c.coeffs[i,10]*y$sumissue + c.coeffs[i,12]*y$black + c.coeffs[i,14]*y$totalmentions + c.coeffs[i,16]*y$pidstrength2  + c.coeffs[i,18]*y$changerdpi + c.coeffs[i,20]*y$reelectcamp)
  
  
  
  aage[y$year[i]>=2000]<-aage+ma-mi
  
  
  
  yed80<-y
  yed80$female[y$year >1980]<-.5581034
  mi<-min(c.coeffs[i,1]+c.coeffs[i,4]*y$education + c.coeffs[i,6]*y$age + c.coeffs[i,8]*yed80$female + c.coeffs[i,10]*y$sumissue + c.coeffs[i,12]*y$black + c.coeffs[i,14]*y$totalmentions + c.coeffs[i,16]*y$pidstrength2  + c.coeffs[i,18]*y$changerdpi + c.coeffs[i,20]*y$reelectcamp)
  yed80$female[y$year >1980]<-.5414663
  ma<-max(c.coeffs[i,1]+c.coeffs[i,4]*y$education + c.coeffs[i,6]*y$age + c.coeffs[i,8]*yed80$female + c.coeffs[i,10]*y$sumissue + c.coeffs[i,12]*y$black + c.coeffs[i,14]*y$totalmentions + c.coeffs[i,16]*y$pidstrength2  + c.coeffs[i,18]*y$changerdpi + c.coeffs[i,20]*y$reelectcamp)
  
  
  afem[y$year[i]>=2000]<-afem+ma-mi
  
  
  
  yed80<-y
  yed80$sumissue[y$year >1980]<-.4429074
  mi<-min(c.coeffs[i,1]+c.coeffs[i,4]*y$education + c.coeffs[i,6]*y$age + c.coeffs[i,8]*y$female + c.coeffs[i,10]*yed80$sumissue + c.coeffs[i,12]*y$black + c.coeffs[i,14]*y$totalmentions + c.coeffs[i,16]*y$pidstrength2  + c.coeffs[i,18]*y$changerdpi + c.coeffs[i,20]*y$reelectcamp)
  yed80$sumissue[y$year >1980]<-.6405873
  ma<-max(c.coeffs[i,1]+c.coeffs[i,4]*y$education + c.coeffs[i,6]*y$age + c.coeffs[i,8]*y$female + c.coeffs[i,10]*yed80$sumissue + c.coeffs[i,12]*y$black + c.coeffs[i,14]*y$totalmentions + c.coeffs[i,16]*y$pidstrength2  + c.coeffs[i,18]*y$changerdpi + c.coeffs[i,20]*y$reelectcamp)
  
  
  asumissue[y$year[i]>=2000]<-asumissue+ma-mi
  
  
  
  yed80<-y
  yed80$black[y$year >1980]<-.5581034
  mi<-min(c.coeffs[i,1]+c.coeffs[i,4]*y$education + c.coeffs[i,6]*y$age + c.coeffs[i,8]*y$female + c.coeffs[i,10]*y$sumissue + c.coeffs[i,12]*yed80$black + c.coeffs[i,14]*y$totalmentions + c.coeffs[i,16]*y$pidstrength2  + c.coeffs[i,18]*y$changerdpi + c.coeffs[i,20]*y$reelectcamp)
  yed80$black[y$year >1980]<-.5414663
  ma<-max(c.coeffs[i,1]+c.coeffs[i,4]*y$education + c.coeffs[i,6]*y$age + c.coeffs[i,8]*y$female + c.coeffs[i,10]*y$sumissue + c.coeffs[i,12]*yed80$black + c.coeffs[i,14]*y$totalmentions + c.coeffs[i,16]*y$pidstrength2  + c.coeffs[i,18]*y$changerdpi + c.coeffs[i,20]*y$reelectcamp)
  
  
  ablack[y$year[i]>=2000]<-ablack+ma-mi
  
  
  
  
  
  yed80<-y
  yed80$totalmentions[y$year >1980]<-8.23119
  mi<-min(c.coeffs[i,1]+c.coeffs[i,4]*y$education + c.coeffs[i,6]*y$age + c.coeffs[i,8]*y$female + c.coeffs[i,10]*y$sumissue + c.coeffs[i,12]*y$black + c.coeffs[i,14]*yed80$totalmentions + c.coeffs[i,16]*y$pidstrength2  + c.coeffs[i,18]*y$changerdpi + c.coeffs[i,20]*y$reelectcamp)
  yed80$totalmentions[y$year >1980]<-8.249172
  ma<-max(c.coeffs[i,1]+c.coeffs[i,4]*y$education + c.coeffs[i,6]*y$age + c.coeffs[i,8]*y$female + c.coeffs[i,10]*y$sumissue + c.coeffs[i,12]*y$black + c.coeffs[i,14]*yed80$totalmentions + c.coeffs[i,16]*y$pidstrength2  + c.coeffs[i,18]*y$changerdpi + c.coeffs[i,20]*y$reelectcamp)
  
  
  atotalmentions[y$year[i]>=2000]<-atotalmentions+ma-mi
  
  
  
  
  yed80<-y
  yed80$pidstrength2[y$year >1980]<-1.198608
  mi<-min(c.coeffs[i,1]+c.coeffs[i,4]*y$education + c.coeffs[i,6]*y$age + c.coeffs[i,8]*y$female + c.coeffs[i,10]*y$sumissue + c.coeffs[i,12]*y$black + c.coeffs[i,14]*yed80$totalmentions + c.coeffs[i,16]*yed80$pidstrength2  + c.coeffs[i,18]*y$changerdpi + c.coeffs[i,20]*y$reelectcamp)
  yed80$pidstrength2[y$year >1980]<-1.204974
  ma<-max(c.coeffs[i,1]+c.coeffs[i,4]*y$education + c.coeffs[i,6]*y$age + c.coeffs[i,8]*y$female + c.coeffs[i,10]*y$sumissue + c.coeffs[i,12]*y$black + c.coeffs[i,14]*y$totalmentions + c.coeffs[i,16]*yed80$pidstrength2  + c.coeffs[i,18]*y$changerdpi + c.coeffs[i,20]*y$reelectcamp)
  
  apidstrength2[y$year[i]>=2000]<-apidstrength2+ma-mi
  
  
  yed80<-y
  yed80$changerdpi[y$year >1980]<-.4429074
  mi<-min(c.coeffs[i,1]+c.coeffs[i,4]*y$education + c.coeffs[i,6]*y$age + c.coeffs[i,8]*y$female + c.coeffs[i,10]*y$sumissue + c.coeffs[i,12]*y$black + c.coeffs[i,14]*yed80$totalmentions + c.coeffs[i,16]*y$pidstrength2  + c.coeffs[i,18]*yed80$changerdpi + c.coeffs[i,20]*y$reelectcamp)
  yed80$changerdpi[y$year >1980]<-.6405873
  ma<-max(c.coeffs[i,1]+c.coeffs[i,4]*y$education + c.coeffs[i,6]*y$age + c.coeffs[i,8]*y$female + c.coeffs[i,10]*y$sumissue + c.coeffs[i,12]*y$black + c.coeffs[i,14]*y$totalmentions + c.coeffs[i,16]*y$pidstrength2  + c.coeffs[i,18]*yed80$changerdpi + c.coeffs[i,20]*y$reelectcamp)
  
  achangerdpi[y$year[i]>=2000]<-achangerdpi+ma-mi
  
  
  yed80<-y
  yed80$reelectcamp[y$year >1980]<-.2857143
  mi<-min(c.coeffs[i,1]+c.coeffs[i,4]*y$education + c.coeffs[i,6]*y$age + c.coeffs[i,8]*y$female + c.coeffs[i,10]*y$sumissue + c.coeffs[i,12]*y$black + c.coeffs[i,14]*yed80$totalmentions + c.coeffs[i,16]*y$pidstrength2  + c.coeffs[i,18]*y$changerdpi + c.coeffs[i,20]*yed80$reelectcamp)
  yed80$reelectcamp[y$year >1980]<-1/3
  ma<-max(c.coeffs[i,1]+c.coeffs[i,4]*y$education + c.coeffs[i,6]*y$age + c.coeffs[i,8]*y$female + c.coeffs[i,10]*y$sumissue + c.coeffs[i,12]*y$black + c.coeffs[i,14]*y$totalmentions + c.coeffs[i,16]*y$pidstrength2  + c.coeffs[i,18]*y$changerdpi + c.coeffs[i,20]*yed80$reelectcamp)
  
  areelectcamp[y$year[i]>=2000]<-areelectcamp+ma-mi
  
  
  print(i)
}

header.names<-c("Ambivalent","Undecided")
var.names<-c("Constant", "Education", "Age", "Female", "Sum Issue", "Black", "Total Mentions", "Pid Strength2", "Changer dpi", "Reelection Campaign")
ambivelence.coeff<-bc.coeff[c(1,4,6,8,10,12,14,16,18,20)]
undecided.coeff<-bc.coeff[c(2,5,7,9,11,13,15,17,19,21)]
ambivelence.se<-bc.se[c(1,4,6,8,10,12,14,16,18,20)]
undecided.se<-bc.se[c(2,5,7,9,11,13,15,17,19,21)]
n<-19639
Elections<-13
ll<- -16046.58
rho<-0.1072
install.packages("xtable")
library(xtable)
xtable(cbind(var.names,round(ambivelence.coeff,6),round(undecided.coeff,6)))



######################################################
######################################################

##Replication Code for Table 3.

######################################################
##Preliminary Steps:

setwd(dirname(file.choose()))

require(foreign)

statadata <- read.dta("modifiedcdf.dta")

head(statadata)

#####################################################
##MNL Estimates Column of Table 3:

#Replacing weights with post-election weights where they exist, per Smidt's .do file.

for(i in 1:nrow(statadata)){
  statadata$weight[i] <- if(is.na(statadata$postweight[i])){statadata$weight[i]}else{statadata$postweight[i]}  
}

summary(statadata$weight)

#Making the dataframe for the mlogit model:

require(mlogit)

#Originally,  0=Repeat Non-voter; 1=Standpatter; 2=Floating Voter; 3=Surger & Decliner.

statadata3 <- statadata[  , c(1, 3, 6, 11, 12, 13, 14, 16, 20, 22, 23, 25, 37, 38, 45)] #this will eliminate unnecessary variables.

statadata3$Newcatvoter2 <- statadata3$catvoter2
statadata3$Newcatvoter2[statadata3$Newcatvoter2== 0] <- 88
statadata3$Newcatvoter2[statadata3$Newcatvoter2== 1] <- 0
statadata3$Newcatvoter2[statadata3$Newcatvoter2== 88] <- 1
table(statadata3$Newcatvoter2) #Check: New frequencies (recoded for 0 and 1 flipped).

frame2 <- mlogit.data(statadata3, shape="wide", choice="Newcatvoter2")

#Modeling Smidt's intermediate mlogit model:

model2 <- mlogit(Newcatvoter2 ~ 0 | efficacy + education + age + female + black + totalmentions + pidstrength2 + sumissue + undecided + ambivalent + changerdpi + reelectcamp, data = frame2, weights = weight)
model2$coefficients

#Check: Coefficients match those in Stata.

#Reweighting process using only data included in model:

head(statadata3)

samp <- na.omit(statadata3)
nrow(samp)

sum <- sum(samp$weight)
samp$weight2 <- samp$weight*(nrow(samp)/sum)

head(samp$weight2)
sum(samp$weight2) 

head(samp)
#save(samp, file = "SampleData.Rdata")
load("SampleData.Rdata")

#Check: These are the same weights as those in Smidt's dataset.

##Another non-bootstrapped model with new weighting scheme:
require(mlogit)
frame3 <- mlogit.data(samp, shape="wide", choice="Newcatvoter2")
model3 <- mlogit(Newcatvoter2 ~ 0 | efficacy + education + age + female + black + totalmentions + pidstrength2 + sumissue + undecided + ambivalent + changerdpi + reelectcamp, data = frame3, weights = weight2)

summary(model3)
model3$coefficients
#1:(intercept)   2:(intercept)   3:(intercept)      1:efficacy      2:efficacy      3:efficacy     1:education     2:education     3:education 
#4.708063498     0.109153862     2.689592543    -1.201483344    -0.057277290    -0.608656858    -0.718535558    -0.161466275    -0.329240137 
#1:age           2:age           3:age        1:female        2:female        3:female         1:black         2:black         3:black 
#-0.047819089    -0.007840548    -0.044745168     0.032702081     0.078483201     0.079334990     0.531217552    -0.342985355     0.354122533 
#1:totalmentions 2:totalmentions 3:totalmentions  1:pidstrength2  2:pidstrength2  3:pidstrength2      1:sumissue      2:sumissue      3:sumissue 
#-0.124261922    -0.005306091    -0.049734345    -0.918418040    -0.741113896    -0.425226354    -0.624280621    -0.423758126    -0.301027260 
#1:undecided     2:undecided     3:undecided    1:ambivalent    2:ambivalent    3:ambivalent    1:changerdpi    2:changerdpi    3:changerdpi 
#0.236252482     0.577622484     0.265172770     0.101892998     0.509656551     0.200991366     0.012324732    -0.013838115     0.015876690 
#1:reelectcamp   2:reelectcamp   3:reelectcamp 
#0.073651851    -0.116921838     0.022827697 

#Check: These coefficients match Smidt's intermediate coefficients.

###############################################################
#Bootstrapping by clusters for MNL coefficients:

#Use "samp" as data.
#use "weight2" as weights.

bootdat <- NULL
onecluster <- NULL
indexdat <- NULL
reps <- 2000
coeffs  <- matrix(NA, nrow = reps, ncol = 39)

set.seed(1004)

clusterboot <- for(i in 1:reps){
  
  bootdat <- NULL  
  model.sim <- NULL
  
  clusters <- names(table(samp$year))
  draws <- sample( 1:length(clusters), length(clusters), replace=TRUE)  
  inds <- clusters[draws]
  print(i)
  
  for(j in 1:length(inds)){
    onecluster <- samp[as.character(samp$year) == inds[j], ]
    
    bootdat <- rbind(bootdat, onecluster)
  }
  
  frame.sim <- mlogit.data(bootdat, choice = "Newcatvoter2", shape = "wide")
  
  model.sim <- tryCatch(mlogit(Newcatvoter2 ~ 0 | efficacy + education + age + female + black + totalmentions + pidstrength2 + sumissue + undecided + ambivalent + changerdpi + reelectcamp, data = frame.sim, weights = weight2),
                        error = function(e){NA})
  
  if (is.na(model.sim)){i = i+1 }else{coeffs.long <- as.data.frame(model.sim$coefficients)[ , 1]}
  
  coeffs[i, ] <- t(as.matrix(coeffs.long[1:39]))
  
}

head(coeffs)
nrow(coeffs)
summary(coeffs)

#save(coeffs, file="BootstrappedCoeffs_Full.Rdata") 

##Attempt to run with more coefficients. ***

####################################################
##Mean of bootstrapped coefficients:

#To begin at this step, 
#load("BootstrappedCoeffs_Full.Rdata")
#head(coeffs)
#dim(coeffs)

coeffs2 <- as.matrix(coeffs)

coeffs2[1771, ] #This is the row that got skipped by the tryCatch() due to colinearity. 
               #The same occurred for Smidt with one of his draws, resulting in 1,999 total draws. 

coeffs2 <- coeffs2[-1771, ] #Eliminating that row from the overall matrix.
head(coeffs2)
dim(coeffs2)

expected.coeffs <- NULL

for(i in 1:ncol(coeffs2)){ 
  vector <-coeffs2[, i]  
  expected.coeffs[i] <- mean(vector)
}

head(expected.coeffs)
expected.coeffs <- as.data.frame(expected.coeffs)
expected.coeffs <- cbind(names(model2$coefficients), expected.coeffs)

####################################################
##Bootstrapped clustered standard errors:

sterrs <- NULL

for(i in 1:ncol(coeffs2)){ 
  vector <-coeffs2[, i]  
  sterrs[i] <- sd(vector)
}

head(sterrs)
sterrs <- as.data.frame(sterrs)
sterrs <- cbind(names(model2$coefficients), sterrs)

#Table of un-bias corrected coefficients and bootstrapped SEs where the 0-2 difference is the Standpatter to Floating Voter MNL estimate.

table <- cbind( names(model2$coefficients), expected.coeffs[ , 2], sterrs[ , 2])
colnames(table) <- c( "" , "expected coefficient", "standard error" )

#save(table, file = "Uncorrected_Coeffs_SEs.Rdata")
#load("Uncorrected_Coeffs_SEs.Rdata")

#Coefficients are very close to those generated by Smidt, but won't be exact due to
#different random number generators across softwares. 

#######################################################################
##calculated p-values by hand.

dim(expected.coeffs)
expected.coeffs <- expected.coeffs[ , -1]

dim(sterrs)
sterrs <- sterrs[ , -1]

##p-values:
for(i in 1:39){
ate <- expected.coeffs[i]
ate.SE <-  sterrs[i]
tstat <- (ate-0)/ate.SE
pvalue[i] <- 2*(1 - pnorm(abs(tstat), 0, 1))
}

pvalue
names(pvalue) <- names(model3$coefficients)


################################################################
## Bias calculations for bootstrapped sample:

#As it is not evident how stata caluclates its bias estimates, I considered bias to be the difference between the
#expected bootstrapped coefficients and the coefficients of the original model. 

bias <- NULL
original.coeffs <- as.data.frame(model2$coefficients)
original.coeffs <- original.coeffs[ , 1]
original.coeffs <- as.vector(original.coeffs[1:39])
expected.vector <- as.vector(expected.coeffs)

bias <- original.coeffs - expected.vector
bias.dat <- as.data.frame(bias)
rownames(bias.dat) <- names(model2$coefficients)
bias.dat

##Vector of bias-coerrected only (thetahat - thetabar):
bc <- expected.vector + as.vector(bias)
names(bc) <- names(model2$coefficients)
bc.SF <- bc[c(2, 5, 8, 11, 14, 17, 20, 23, 26, 29, 32, 35, 38)]
as.data.frame(bc.SF)

##After getting this bias-correction statistic for each coefficient, I apply it to all bootstrapped coefficients and take their mean. 
##This allows the use of a bias-corrected set of 2,000 rows of bootstrapped coefficients for use in calculating across-era effects. 

#use coeffs2
#use bias.dat

head(coeffs2)
dim(coeffs2)
bc.coeffs2 <- coeffs2
head(bc.coeffs2)
dim(bc.coeffs2)

bias.dat2 <- t(bias.dat)
head(bias.dat2)
dim(bias.dat2)

for(i in 1:1999){
  bc.coeffs2[i, ] <-  bc.coeffs2[i, ] + bias.dat2
}

head(bc.coeffs2)

##Getting means of these bc coefficients for use in Table 3:
bc.expected.coeffs <- NULL

for(i in 1:ncol(bc.coeffs2)){ 
  vector <-bc.coeffs2[, i]  
  bc.expected.coeffs[i] <- mean(vector)
}

#Extracting only those for Standpatter-Floating:
bcExp <- bc.expected.coeffs[ c(2, 5, 8, 11, 14, 17, 20, 23, 26, 29, 32, 35, 38)]
names(bcExp)<- c("intercept", "efficacy", "education", "age", "female", "back", "totalmentions", "pidstrength2", 
                 "sumissue", "undecided", "ambivalent", "changerpdi", "reelectcamp")

#save(bc.expected.coeffs, file = "BC_Coeffs_All_32016.Rdata")  
#To begin from this point, load "BC_Coeffs_All_32016.Rdata".

###############################################################################################
###############################################################################################
########Calculating Across-Era Effects (columns 2 & 3 of Table 3):

##Preparing original sample values and creating backup variables:

part.samp <- samp[ , c(1, 9, 4, 5, 6, 7, 13, 8, 15, 3, 14, 12, 11, 16, 17)] 
intercept <- rep(1, 14424) 
part.samp <- cbind(intercept, part.samp) 

names(part.samp) <- c("intercept", "year", "efficacy", "education", "age", "female", "black", "totalmentions", "pidstrength2", 
                       "sumissue", "undecided", "ambivalent", "changerpdi", "reelectcamp", "Newcatvoter2", "weight2")
dim(part.samp)                        
head(part.samp) 

#part.samp$bkefficacy <- part.samp$efficacy  ##Creating backup variables to backfill means after each use. 
#part.samp$bkeducation <- part.samp$education
#part.samp$bkage <- part.samp$age
#part.samp$bkfemale <- part.samp$female
#part.samp$bkblack <- part.samp$black
#part.samp$bktotalmentions <- part.samp$totalmentions
#part.samp$bkpidstrength2 <- part.samp$pidstrength2
#part.samp$bksumissue <- part.samp$sumissue
#part.samp$bkundecided <- part.samp$undecided
#part.samp$bkambivalent <- part.samp$ambivalent
#part.samp$bkchangerpdi <- part.samp$changerpdi
#part.samp$bkreelectcamp <- part.samp$reelectcamp

#################################################
##Relevant Data for use from this point ( skip these steps if data is already loaded. 

head(part.samp) #sample with relevant variables only.
dim(part.samp)
#save(part.samp, file = "SampleData_Clean.Rdata")
#To begin from this point, use:
load("SampleData_Clean.Rdata") # for sample.


#To begin from this point, use:
load("BootstrappedCoeffs_Full.Rdata")  #for coefficients.
#Then,
coeffs2 <- as.matrix(coeffs)
#coeffs2[1771, ]   #Skip this step if using coeffs2 from above.
#coeffs2 <- coeffs2[-1771, ] #Eliminating that row from the overall matrix.
head(coeffs2)
dim(coeffs2)

#####################################
##creating new variables for storage (this step can be skipped if user prefers to create a new object each time):

temp <- matrix(NA, nrow=14424, ncol = 64)
temp <- as.data.frame(temp)
head(temp)
names(temp) <- c("d_educ0", "d_educ1", "d_educ2", "d_educ3", 
                 "d_age0", "d_age1", "d_age2", "d_age3",
                 "d_fem0", "d_fem1", "d_fem2", "d_fem3",
                 "d_black0", "d_black1", "d_black2", "d_black3",
                 "d_demo0", "d_demo1", "d_demo2", "d_demo3",
                 "d_eff0", "d_eff1", "d_eff2", "d_eff3",
                 "d_tm0", "d_tm1", "d_tm2", "d_tm3",
                 "d_strength0", "d_strength1", "d_strength2", "d_strength3",
                 "d_engage0", "d_engage1", "d_engage2", "d_engage3",
                 "d_change0", "d_change1", "d_change2", "d_change3",
                 "d_reelect0", "d_reelect1","d_reelect2", "d_reelect3",
                 "d_election0", "d_election1", "d_election2", "d_election3",
                 "d_issue0", "d_issue1", "d_issue2", "d_issue3",
                 "d_und0", "d_und1", "d_und2", "d_und3",
                 "d_amb0", "d_amb1", "d_amb2", "d_amb3",
                 "d_clarity0", "d_clarity1", "d_clarity2", "d_clarity3")

## Appending new variables to part.samp data:

combo.dat <- cbind(part.samp, temp)
dim(combo.dat)
head(combo.dat)

##Getting betas matrix formatted for use in across-era effects:

bcsim.betas1 <- bc.coeffs2[ , c(1, 4, 7, 10, 13, 16, 19, 22, 25, 28, 31, 34, 37)]
dim(bcsim.betas1)
bcsim.betas2 <- bc.coeffs2[ , c(2, 5, 8, 11, 14, 17, 20, 23, 26, 29, 32, 35, 38)]
bcsim.betas3 <- bc.coeffs2[ , c(3, 6, 9, 12, 15, 18, 21, 24, 27, 30, 33, 36, 39)]
bcsim.betas.base <- matrix(data = 0, nrow = 1999, ncol = 13)
colnames(bcsim.betas.base) <- c("0:intercept", "0:efficacy", "0:education", "0:age", "0:female", "0:back", "0:totalmentions", "0:pidstrength2", 
                              "0:sumissue", "0:undecided", "0:ambivalent", "0:changerpdi", "0:reelectcamp")
bcsimbetas <- cbind( bcsim.betas1, bcsim.betas.base , bcsim.betas2, bcsim.betas3)  ##Putting base of 0s in to match Smidt's bc coefficient matrix.
bcsimbetas <- as.matrix(bcsimbetas) 
dim(bcsimbetas)

##########################################

##For across-era effects, hold each variable constant at weighted mean pre-1980 and weighted mean since 2000 in turn. 
##Keep other variables at observed values.
##Calculate differences across eras in predicted probability of being floating voter (versus baseline of standpatter (0 in my data)).

##Caluclating weighted means:

nrow(combo.dat)

wmeans.80 <- rep(NA, 12)

wmeans.80[1] <- weighted.mean(combo.dat$efficacy[combo.dat$year < 1980], combo.dat$weight2[combo.dat$year < 1980])
wmeans.80[2] <- weighted.mean(combo.dat$education[combo.dat$year < 1980],  combo.dat$weight2[combo.dat$year < 1980])
wmeans.80[3] <- weighted.mean(combo.dat$age[combo.dat$year < 1980],  combo.dat$weight2[combo.dat$year < 1980])
wmeans.80[4] <- weighted.mean(combo.dat$female[combo.dat$year < 1980],  combo.dat$weight2[combo.dat$year < 1980])
wmeans.80[5] <- weighted.mean(combo.dat$black[combo.dat$year < 1980],  combo.dat$weight2[combo.dat$year < 1980])
wmeans.80[6] <- weighted.mean(combo.dat$totalmentions[combo.dat$year < 1980],  combo.dat$weight2[combo.dat$year < 1980])
wmeans.80[7] <- weighted.mean(combo.dat$pidstrength2[combo.dat$year < 1980],  combo.dat$weight2[combo.dat$year < 1980])
wmeans.80[8] <- weighted.mean(combo.dat$sumissue[combo.dat$year < 1980], combo.dat$weight2[combo.dat$year < 1980])
wmeans.80[9] <- weighted.mean(combo.dat$undecided[combo.dat$year < 1980], combo.dat$weight2[combo.dat$year < 1980])
wmeans.80[10] <- weighted.mean(combo.dat$ambivalent[combo.dat$year < 1980],  combo.dat$weight2[combo.dat$year < 1980])
wmeans.80[11] <- weighted.mean(combo.dat$changerpdi[combo.dat$year < 1980],  combo.dat$weight2[combo.dat$year < 1980])
wmeans.80[12] <- weighted.mean(combo.dat$reelectcamp[combo.dat$year < 1980],  combo.dat$weight2[combo.dat$year < 1980])

sum(combo.dat$weight2[combo.dat$year < 1980])

wmeans.80[8] #0.443376  (<1980)
##Compare to Smidt mean for pre-1980 era, .4429074. 

wmeans.80
# 0.64464028  1.16874273 45.16468369  0.56588710  0.09259069  8.47937489  1.20975875  0.44337601  0.11687445  0.26519622  3.00107789  0.33315155
wmeans.80 <- c(1, wmeans.80)
names(wmeans.80) <- c("intercept", "efficacy", "education", "age", "female", "black", "totalmentions", "pidstrength2", 
                      "sumissue", "undecided", "ambivalent", "changerpdi", "reelectcamp")

##Getting post-2000 weighted means for all variables:
##By the same logic as above, I will use the na.omit sample only, with weights. 

wmeans.00 <- rep(NA, 12)

wmeans.00[1] <- weighted.mean(combo.dat$efficacy[combo.dat$year > 1999], combo.dat$weight2[combo.dat$year > 1999])
wmeans.00[2] <- weighted.mean(combo.dat$education[combo.dat$year > 1999], combo.dat$weight2[combo.dat$year > 1999])  
wmeans.00[3] <- weighted.mean(combo.dat$age[combo.dat$year > 1999], combo.dat$weight2[combo.dat$year > 1999])
wmeans.00[4] <- weighted.mean(combo.dat$female[combo.dat$year > 1999], combo.dat$weight2[combo.dat$year > 1999])
wmeans.00[5] <- weighted.mean(combo.dat$black[combo.dat$year > 1999], combo.dat$weight2[combo.dat$year > 1999])
wmeans.00[6] <- weighted.mean(combo.dat$totalmentions[combo.dat$year > 1999], combo.dat$weight2[combo.dat$year > 1999])
wmeans.00[7] <- weighted.mean(combo.dat$pidstrength2[combo.dat$year > 1999], combo.dat$weight2[combo.dat$year > 1999])
wmeans.00[8] <- weighted.mean(combo.dat$sumissue[combo.dat$year > 1999], combo.dat$weight2[combo.dat$year > 1999])
wmeans.00[9] <- weighted.mean(combo.dat$undecided[combo.dat$year > 1999], combo.dat$weight2[combo.dat$year > 1999])
wmeans.00[10] <- weighted.mean(combo.dat$ambivalent[combo.dat$year > 1999], combo.dat$weight2[combo.dat$year > 1999])
wmeans.00[11] <- weighted.mean(combo.dat$changerpdi[combo.dat$year > 1999], combo.dat$weight2[combo.dat$year > 1999])
wmeans.00[12] <- weighted.mean(combo.dat$reelectcamp[combo.dat$year > 1999], combo.dat$weight2[combo.dat$year > 1999])

wmeans.00
#  0.47515012  1.75089365 46.97885995  0.54402266  0.15855633  8.55806766  1.25158748  0.62636900  0.04487127  0.18708400  3.40331192  0.42811759
#Compare mean.00[8]: 0.626369 to Smidt's post-2000 sumissue mean of .6405873. Ok, and a bit better than nonweighted. 

wmeans.00 <- c(1, wmeans.00)
names(wmeans.00) <- c("intercept", "efficacy", "education", "age", "female", "black", "totalmentions", "pidstrength2", 
                      "sumissue", "undecided", "ambivalent", "changerpdi", "reelectcamp")

##While not exact, these weighted means using only the sample are closer to Smidt's means
##than other specifications. See Alternative Specifications script for non-weighted means, means from full dataset, etc. . 

###########################################
##Effects calculations for each variable (single):

##########
##Sumisse:
calcdata80 <- (combo.dat[ , c(1, 3:14)])
wmeans.80[9]  #using weighted mean for the variable in question.
calcdata80[ , 9] <- rep(wmeans.80[9], nrow(calcdata80))
calcdata80 <- as.matrix(calcdata80)
head(calcdata80)

calcdata00 <- (combo.dat[ , c(1, 3:14)])
wmeans.00[9]
calcdata00[ , 9] <- rep(wmeans.00[9], nrow(calcdata00))
calcdata00 <- as.matrix(calcdata00)
head(calcdata00)

d_sumissue2 <- NULL  #We only really need to recreate Standpatter to Floating Voters effect:

for(i in 1:1999){
  effect_sumissue.80.SP_FL <- exp(calcdata80 %*% bcsimbetas[i, 27:39])/(1 + exp(calcdata80 %*% bcsimbetas[i, 1:13]) + exp(calcdata80%*% bcsimbetas[i, 27:39]) +  exp(calcdata80 %*% bcsimbetas[i, 40:52]))
  effect_sumissue.00.SP_FL <- exp(calcdata00 %*% bcsimbetas[i, 27:39])/(1 + exp(calcdata00 %*% bcsimbetas[i, 1:13]) + exp(calcdata00 %*% bcsimbetas[i, 27:39]) +  exp(calcdata00 %*% bcsimbetas[i, 40:52]))
  diff <- effect_sumissue.00.SP_FL - effect_sumissue.80.SP_FL    
  d_sumissue2[i] <- mean(diff)
}

head(d_sumissue2)
length(d_sumissue2)
mean(d_sumissue2)
#[1] -0.00433171  #weighted means tend to underestimate, whereas unweighted means overestimate. Use weighted means and bc sims. 
sd(d_sumissue2)
#[1]  0.001845083

##p-values:
ate <- mean(d_sumissue2)
ate.SE <-  sd(d_sumissue2)
tstat <- (ate-0)/ate.SE
pvalue <- 2*(1 - pnorm(abs(tstat), 0, 1))
# 0.01888949


#############
##Ambivalent:

x <- 11

calcdata80 <- (combo.dat[ , c(1, 3:14)])
wmeans.80[x]
calcdata80[ , x] <- rep(wmeans.80[x], nrow(calcdata80))
calcdata80 <- as.matrix(calcdata80)
head(calcdata80)

calcdata00 <- (combo.dat[ , c(1, 3:14)])
wmeans.00[x]
calcdata00[ , x ] <- rep(wmeans.00[x], nrow(calcdata00))
calcdata00 <- as.matrix(calcdata00)
head(calcdata00)

d_ambivalent2 <- NULL

for(i in 1:1999){
  effect_undecided.80.SP_FL <- exp(calcdata80 %*% bcsimbetas[i, 27:39])/(1 + exp(calcdata80 %*% bcsimbetas[i, 1:13]) + exp(calcdata80%*% bcsimbetas[i, 27:39]) +  exp(calcdata80 %*% bcsimbetas[i, 40:52]))
  effect_undecided.00.SP_FL <- exp(calcdata00 %*% bcsimbetas[i, 27:39])/(1 + exp(calcdata00 %*% bcsimbetas[i, 1:13]) + exp(calcdata00 %*% bcsimbetas[i, 27:39]) +  exp(calcdata00 %*% bcsimbetas[i, 40:52]))
  diff <- effect_undecided.00.SP_FL - effect_undecided.80.SP_FL    
  d_ambivalent2[i] <- mean(diff)
}

head(d_ambivalent2)
length(d_ambivalent2)
mean(d_ambivalent2)
# -0.003093529 
sd(d_ambivalent2)
# 0.0004560158

##p-values:
ate <- mean(d_ambivalent2)
ate.SE <-  sd(d_ambivalent2)
tstat <- (ate-0)/ate.SE
pvalue <- 2*(1 - pnorm(abs(tstat), 0, 1))
#1.170397e-11


###########
##Undecided:

calcdata80 <- (combo.dat[ , c(1, 3:14)])
wmeans.80[10]
calcdata80[ ,10] <- rep(wmeans.80[10], nrow(calcdata80))
calcdata80 <- as.matrix(calcdata80)
head(calcdata80)

calcdata00 <- (combo.dat[ , c(1, 3:14)])
wmeans.00[10]
calcdata00[ , 10] <- rep(wmeans.00[10], nrow(calcdata00))
calcdata00 <- as.matrix(calcdata00)
head(calcdata00)

d_undecided2 <- NULL

for(i in 1:1999){
  effect_undecided.80.SP_FL <- exp(calcdata80 %*% bcsimbetas[i, 27:39])/(1 + exp(calcdata80 %*% bcsimbetas[i, 1:13]) + exp(calcdata80%*% bcsimbetas[i, 27:39]) +  exp(calcdata80 %*% bcsimbetas[i, 40:52]))
  effect_undecided.00.SP_FL <- exp(calcdata00 %*% bcsimbetas[i, 27:39])/(1 + exp(calcdata00 %*% bcsimbetas[i, 1:13]) + exp(calcdata00 %*% bcsimbetas[i, 27:39]) +  exp(calcdata00 %*% bcsimbetas[i, 40:52]))
  diff <- effect_undecided.00.SP_FL - effect_undecided.80.SP_FL    
  d_undecided2[i] <- mean(diff)
}

head(d_undecided2)
length(d_undecided2)
mean(d_undecided2)
#  -0.003078924
sd(d_undecided2)
# 0.0005288059

##p-values:
ate <- mean(d_undecided2)
ate.SE <-  sd(d_undecided2)
tstat <- (ate-0)/ate.SE
pvalue <- 2*(1 - pnorm(abs(tstat), 0, 1))
#5.800593e-09

################
##Efficacy:
x <- 2

calcdata80 <- (combo.dat[ , c(1, 3:14)])
wmeans.80[x]
calcdata80[ , x] <- rep(wmeans.80[x], nrow(calcdata80))
calcdata80 <- as.matrix(calcdata80)
head(calcdata80)

calcdata00 <- (combo.dat[ , c(1, 3:14)])
wmeans.00[x]
calcdata00[ , x ] <- rep(wmeans.00[x], nrow(calcdata00))
calcdata00 <- as.matrix(calcdata00)
head(calcdata00)

d_efficacy2 <- NULL

for(i in 1:1999){
  effect_undecided.80.SP_FL <- exp(calcdata80 %*% bcsimbetas[i, 27:39])/(1 + exp(calcdata80 %*% bcsimbetas[i, 1:13]) + exp(calcdata80%*% bcsimbetas[i, 27:39]) +  exp(calcdata80 %*% bcsimbetas[i, 40:52]))
  effect_undecided.00.SP_FL <- exp(calcdata00 %*% bcsimbetas[i, 27:39])/(1 + exp(calcdata00 %*% bcsimbetas[i, 1:13]) + exp(calcdata00 %*% bcsimbetas[i, 27:39]) +  exp(calcdata00 %*% bcsimbetas[i, 40:52]))
  diff <- effect_undecided.00.SP_FL - effect_undecided.80.SP_FL    
  d_efficacy2[i] <- mean(diff)
}

head(d_efficacy2)
length(d_efficacy2)
mean(d_efficacy2)
# -0.004727966
sd(d_efficacy2)
# 0.001601033

##p-values:
ate <- mean(d_efficacy2)
ate.SE <-  sd(d_efficacy2)
tstat <- (ate-0)/ate.SE
pvalue <- 2*(1 - pnorm(abs(tstat), 0, 1))
#0.003146286


################
##Education:
x <- 3

calcdata80 <- (combo.dat[ , c(1, 3:14)])
wmeans.80[x]
calcdata80[ , x] <- rep(wmeans.80[x], nrow(calcdata80))
calcdata80 <- as.matrix(calcdata80)
head(calcdata80)

calcdata00 <- (combo.dat[ , c(1, 3:14)])
wmeans.00[x]
calcdata00[ , x ] <- rep(wmeans.00[x], nrow(calcdata00))
calcdata00 <- as.matrix(calcdata00)
head(calcdata00)

d_education2 <- NULL

for(i in 1:1999){
  effect_undecided.80.SP_FL <- exp(calcdata80 %*% bcsimbetas[i, 27:39])/(1 + exp(calcdata80 %*% bcsimbetas[i, 1:13]) + exp(calcdata80%*% bcsimbetas[i, 27:39]) +  exp(calcdata80 %*% bcsimbetas[i, 40:52]))
  effect_undecided.00.SP_FL <- exp(calcdata00 %*% bcsimbetas[i, 27:39])/(1 + exp(calcdata00 %*% bcsimbetas[i, 1:13]) + exp(calcdata00 %*% bcsimbetas[i, 27:39]) +  exp(calcdata00 %*% bcsimbetas[i, 40:52]))
  diff <- effect_undecided.00.SP_FL - effect_undecided.80.SP_FL    
  d_education2[i] <- mean(diff)
}

head(d_education2)
length(d_education2)
mean(d_education2)
# 0.002195335
sd(d_education2)
# 0.001538582

##p-values:
ate <- mean(d_education2)
ate.SE <-  sd(d_education2)
tstat <- (ate-0)/ate.SE
pvalue <- 2*(1 - pnorm(abs(tstat), 0, 1))
#0.1536214


################
##Age:
x <- 4

calcdata80 <- (combo.dat[ , c(1, 3:14)])
wmeans.80[x]
calcdata80[ , x] <- rep(wmeans.80[x], nrow(calcdata80))
calcdata80 <- as.matrix(calcdata80)
head(calcdata80)

calcdata00 <- (combo.dat[ , c(1, 3:14)])
wmeans.00[x]
calcdata00[ , x ] <- rep(wmeans.00[x], nrow(calcdata00))
calcdata00 <- as.matrix(calcdata00)
head(calcdata00)

d_age2 <- NULL

for(i in 1:1999){
  effect_undecided.80.SP_FL <- exp(calcdata80 %*% bcsimbetas[i, 27:39])/(1 + exp(calcdata80 %*% bcsimbetas[i, 1:13]) + exp(calcdata80%*% bcsimbetas[i, 27:39]) +  exp(calcdata80 %*% bcsimbetas[i, 40:52]))
  effect_undecided.00.SP_FL <- exp(calcdata00 %*% bcsimbetas[i, 27:39])/(1 + exp(calcdata00 %*% bcsimbetas[i, 1:13]) + exp(calcdata00 %*% bcsimbetas[i, 27:39]) +  exp(calcdata00 %*% bcsimbetas[i, 40:52]))
  diff <- effect_undecided.00.SP_FL - effect_undecided.80.SP_FL    
  d_age2[i] <- mean(diff)
}

head(d_age2)
length(d_age2)
mean(d_age2)   ##Need to revisit this one. Check against Smidt's means.
# 0.001826073
sd(d_age2)
# 0.0003009961

##p-values:
ate <- mean(d_age2)
ate.SE <-  sd(d_age2)
tstat <- (ate-0)/ate.SE
pvalue <- 2*(1 - pnorm(abs(tstat), 0, 1))
#1.305119e-09


################
##Female:
x <- 5

calcdata80 <- (combo.dat[ , c(1, 3:14)])
wmeans.80[x]
calcdata80[ , x] <- rep(wmeans.80[x], nrow(calcdata80))
calcdata80 <- as.matrix(calcdata80)
head(calcdata80)

calcdata00 <- (combo.dat[ , c(1, 3:14)])
wmeans.00[x]
calcdata00[ , x] <- rep(wmeans.00[x], nrow(calcdata00))
calcdata00 <- as.matrix(calcdata00)
head(calcdata00)

d_female2 <- NULL

for(i in 1:1999){
  effect_undecided.80.SP_FL <- exp(calcdata80 %*% bcsimbetas[i, 27:39])/(1 + exp(calcdata80 %*% bcsimbetas[i, 1:13]) + exp(calcdata80%*% bcsimbetas[i, 27:39]) +  exp(calcdata80 %*% bcsimbetas[i, 40:52]))
  effect_undecided.00.SP_FL <- exp(calcdata00 %*% bcsimbetas[i, 27:39])/(1 + exp(calcdata00 %*% bcsimbetas[i, 1:13]) + exp(calcdata00 %*% bcsimbetas[i, 27:39]) +  exp(calcdata00 %*% bcsimbetas[i, 40:52]))
  diff <- effect_undecided.00.SP_FL - effect_undecided.80.SP_FL    
  d_female2[i] <- mean(diff)
}

head(d_female2)
length(d_female2)
mean(d_female2)
# -0.0001109187
sd(d_female2)
#  0.0001150836

##p-values:
ate <- mean(d_female2)
ate.SE <-  sd(d_female2)
tstat <- (ate-0)/ate.SE
pvalue <- 2*(1 - pnorm(abs(tstat), 0, 1))
#0.3351414


################
##Black:
x <- 6

calcdata80 <- (combo.dat[ , c(1, 3:14)])
wmeans.80[x]
calcdata80[ , x] <- rep(wmeans.80[x], nrow(calcdata80))
calcdata80 <- as.matrix(calcdata80)
head(calcdata80)

calcdata00 <- (combo.dat[ , c(1, 3:14)])
wmeans.00[x]
calcdata00[ , x] <- rep(wmeans.00[x], nrow(calcdata00))
calcdata00 <- as.matrix(calcdata00)
head(calcdata00)

d_black2 <- NULL

for(i in 1:1999){
  effect_undecided.80.SP_FL <- exp(calcdata80 %*% bcsimbetas[i, 27:39])/(1 + exp(calcdata80 %*% bcsimbetas[i, 1:13]) + exp(calcdata80%*% bcsimbetas[i, 27:39]) +  exp(calcdata80 %*% bcsimbetas[i, 40:52]))
  effect_undecided.00.SP_FL <- exp(calcdata00 %*% bcsimbetas[i, 27:39])/(1 + exp(calcdata00 %*% bcsimbetas[i, 1:13]) + exp(calcdata00 %*% bcsimbetas[i, 27:39]) +  exp(calcdata00 %*% bcsimbetas[i, 40:52]))
  diff <- effect_undecided.00.SP_FL - effect_undecided.80.SP_FL    
  d_black2[i] <- mean(diff)
}

head( d_black2)
length( d_black2)
mean( d_black2)
#  -0.003196287
sd( d_black2)
#  0.001301696

##p-values:
ate <- mean(d_black2)
ate.SE <-  sd(d_black2)
tstat <- (ate-0)/ate.SE
pvalue <- 2*(1 - pnorm(abs(tstat), 0, 1))
#0.01406973


################
##Total Mentions:
x <- 7

calcdata80 <- (combo.dat[ , c(1, 3:14)])
wmeans.80[x]
calcdata80[ , x] <- rep(wmeans.80[x], nrow(calcdata80))
calcdata80 <- as.matrix(calcdata80)
head(calcdata80)

calcdata00 <- (combo.dat[ , c(1, 3:14)])
wmeans.00[x]
calcdata00[ , x] <- rep(wmeans.00[x], nrow(calcdata00))
calcdata00 <- as.matrix(calcdata00)
head(calcdata00)

d_totalmentions2 <- NULL

for(i in 1:1999){
  effect_undecided.80.SP_FL <- exp(calcdata80 %*% bcsimbetas[i, 27:39])/(1 + exp(calcdata80 %*% bcsimbetas[i, 1:13]) + exp(calcdata80%*% bcsimbetas[i, 27:39]) +  exp(calcdata80 %*% bcsimbetas[i, 40:52]))
  effect_undecided.00.SP_FL <- exp(calcdata00 %*% bcsimbetas[i, 27:39])/(1 + exp(calcdata00 %*% bcsimbetas[i, 1:13]) + exp(calcdata00 %*% bcsimbetas[i, 27:39]) +  exp(calcdata00 %*% bcsimbetas[i, 40:52]))
  diff <- effect_undecided.00.SP_FL - effect_undecided.80.SP_FL    
  d_totalmentions2[i] <- mean(diff)
}

head(d_totalmentions2)
length(d_totalmentions2)
mean(d_totalmentions2)
#  0.0002150665
sd(d_totalmentions2)
#  4.744994e-05

##p-values:
ate <- mean(d_totalmentions2)
ate.SE <-  sd(d_totalmentions2)
tstat <- (ate-0)/ate.SE
pvalue <- 2*(1 - pnorm(abs(tstat), 0, 1))
#5.829209e-06


################
##Pidstrength2
x <- 8

calcdata80 <- (combo.dat[ , c(1, 3:14)])
wmeans.80[x]
calcdata80[ , x] <- rep(wmeans.80[x], nrow(calcdata80))
calcdata80 <- as.matrix(calcdata80)
head(calcdata80)

calcdata00 <- (combo.dat[ , c(1, 3:14)])
wmeans.00[x]
calcdata00[ , x] <- rep(wmeans.00[x], nrow(calcdata00))
calcdata00 <- as.matrix(calcdata00)
head(calcdata00)

d_pidstrength2 <- NULL

for(i in 1:1999){
  effect_undecided.80.SP_FL <- exp(calcdata80 %*% bcsimbetas[i, 27:39])/(1 + exp(calcdata80 %*% bcsimbetas[i, 1:13]) + exp(calcdata80%*% bcsimbetas[i, 27:39]) +  exp(calcdata80 %*% bcsimbetas[i, 40:52]))
  effect_undecided.00.SP_FL <- exp(calcdata00 %*% bcsimbetas[i, 27:39])/(1 + exp(calcdata00 %*% bcsimbetas[i, 1:13]) + exp(calcdata00 %*% bcsimbetas[i, 27:39]) +  exp(calcdata00 %*% bcsimbetas[i, 40:52]))
  diff <- effect_undecided.00.SP_FL - effect_undecided.80.SP_FL    
  d_pidstrength2[i] <- mean(diff)
}

head(d_pidstrength2)
length(d_pidstrength2)
mean(d_pidstrength2)
# -0.001962537
sd(d_pidstrength2)
#  0.0001949428

##p-values:
ate <- mean(d_pidstrength2)
ate.SE <-  sd(d_pidstrength2)
tstat <- (ate-0)/ate.SE
pvalue <- 2*(1 - pnorm(abs(tstat), 0, 1))
# 0


################
##Reelectcamp
x <- 12

calcdata80 <- (combo.dat[ , c(1, 3:14)])
wmeans.80[x]
calcdata80[ , x] <- rep(wmeans.80[x], nrow(calcdata80))
calcdata80 <- as.matrix(calcdata80)
head(calcdata80)

calcdata00 <- (combo.dat[ , c(1, 3:14)])
wmeans.00[x]
calcdata00[ , x] <- rep(wmeans.00[x], nrow(calcdata00))
calcdata00 <- as.matrix(calcdata00)
head(calcdata00)

d_changerpdi2 <- NULL

for(i in 1:1999){
  effect_undecided.80.SP_FL <- exp(calcdata80 %*% bcsimbetas[i, 27:39])/(1 + exp(calcdata80 %*% bcsimbetas[i, 1:13]) + exp(calcdata80%*% bcsimbetas[i, 27:39]) +  exp(calcdata80 %*% bcsimbetas[i, 40:52]))
  effect_undecided.00.SP_FL <- exp(calcdata00 %*% bcsimbetas[i, 27:39])/(1 + exp(calcdata00 %*% bcsimbetas[i, 1:13]) + exp(calcdata00 %*% bcsimbetas[i, 27:39]) +  exp(calcdata00 %*% bcsimbetas[i, 40:52]))
  diff <- effect_undecided.00.SP_FL - effect_undecided.80.SP_FL    
  d_changerpdi2 [i] <- mean(diff)
}

head(d_changerpdi2)
length(d_changerpdi2)
mean(d_changerpdi2)
#  -0.0007492672
sd(d_changerpdi2)
#  0.001884253

##p-values:
ate <- mean(d_changerpdi2)
ate.SE <-  sd(d_changerpdi2)
tstat <- (ate-0)/ate.SE
pvalue <- 2*(1 - pnorm(abs(tstat), 0, 1))
# 0.6908905

################
##Reelectcamp
x <- 13

calcdata80 <- (combo.dat[ , c(1, 3:14)])
wmeans.80[x]
calcdata80[ , x] <- rep(wmeans.80[x], nrow(calcdata80))
calcdata80 <- as.matrix(calcdata80)
head(calcdata80)

calcdata00 <- (combo.dat[ , c(1, 3:14)])
wmeans.00[x]
calcdata00[ , x] <- rep(wmeans.00[x], nrow(calcdata00))
calcdata00 <- as.matrix(calcdata00)
head(calcdata00)

d_reelectcamp2 <- NULL

for(i in 1:1999){
  effect_undecided.80.SP_FL <- exp(calcdata80 %*% bcsimbetas[i, 27:39])/(1 + exp(calcdata80 %*% bcsimbetas[i, 1:13]) + exp(calcdata80%*% bcsimbetas[i, 27:39]) +  exp(calcdata80 %*% bcsimbetas[i, 40:52]))
  effect_undecided.00.SP_FL <- exp(calcdata00 %*% bcsimbetas[i, 27:39])/(1 + exp(calcdata00 %*% bcsimbetas[i, 1:13]) + exp(calcdata00 %*% bcsimbetas[i, 27:39]) +  exp(calcdata00 %*% bcsimbetas[i, 40:52]))
  diff <- effect_undecided.00.SP_FL - effect_undecided.80.SP_FL    
  d_reelectcamp2  [i] <- mean(diff)
}

head(d_reelectcamp2)
length(d_reelectcamp2)
mean(d_reelectcamp2)
#  -0.001288704
sd(d_reelectcamp2)
#  0.001434291

##p-values:
ate <- mean(d_reelectcamp2)
ate.SE <-  sd(d_reelectcamp2)
tstat <- (ate-0)/ate.SE
pvalue <- 2*(1 - pnorm(abs(tstat), 0, 1))
# 0.3689214



#####################################
######################################
##Combined across era effects: 

##For these calculations, Smidt holds all of the variables in a category (see table 3 for categories) at their mean at once, and then 
##conducts the same caluclation of across-era effect.

##########
##Clarity of Choice: (sumissue(9), ambivalent(11), undecided(10) )

calcdata80 <- (combo.dat[ , c(1, 3:14)])
wmeans.80[9]  #using weighted means for the variables in question.
wmeans.80[11]
wmeans.80[10]
calcdata80[ , 9] <- rep(wmeans.80[9], nrow(calcdata80))
calcdata80[ , 11] <- rep(wmeans.80[11], nrow(calcdata80))
calcdata80[ , 10] <- rep(wmeans.80[10], nrow(calcdata80))
calcdata80 <- as.matrix(calcdata80)
head(calcdata80)

calcdata00 <- (combo.dat[ , c(1, 3:14)])
wmeans.00[9]
wmeans.00[11]
wmeans.00[10]
calcdata00[ , 9] <- rep(wmeans.00[9], nrow(calcdata00))
calcdata00[ , 11] <- rep(wmeans.00[11], nrow(calcdata00))
calcdata00[ , 10] <- rep(wmeans.00[10], nrow(calcdata00))
calcdata00 <- as.matrix(calcdata00)
head(calcdata00)

d_clarity2 <- NULL  #We only really need to recreate Standpatter to Floating Voters effect:

for(i in 1:1999){
  effect_sumissue.80.SP_FL <- exp(calcdata80 %*% bcsimbetas[i, 27:39])/(1 + exp(calcdata80 %*% bcsimbetas[i, 1:13]) + exp(calcdata80%*% bcsimbetas[i, 27:39]) +  exp(calcdata80 %*% bcsimbetas[i, 40:52]))
  effect_sumissue.00.SP_FL <- exp(calcdata00 %*% bcsimbetas[i, 27:39])/(1 + exp(calcdata00 %*% bcsimbetas[i, 1:13]) + exp(calcdata00 %*% bcsimbetas[i, 27:39]) +  exp(calcdata00 %*% bcsimbetas[i, 40:52]))
  diff <- effect_sumissue.00.SP_FL - effect_sumissue.80.SP_FL    
  d_clarity2[i] <- mean(diff)
}

head(d_clarity2)
length(d_clarity2)
mean(d_clarity2)
-0.01039399
sd(d_clarity2)
0.001870485

##p-values:
ate <- mean(d_clarity2)
ate.SE <-  sd(d_clarity2)
tstat <- (ate-0)/ate.SE
pvalue <- 2*(1 - pnorm(abs(tstat), 0, 1))
# 2.746977e-08


##########
##General Political Orientation: (efficacy(2), totalmentions(7), pidstrength2(8))

X <- 2
Y <- 7
Z <- 8

calcdata80 <- (combo.dat[ , c(1, 3:14)])
wmeans.80[X]  
wmeans.80[Y]
wmeans.80[Z]
calcdata80[ , X] <- rep(wmeans.80[X], nrow(calcdata80))
calcdata80[ , Y] <- rep(wmeans.80[Y], nrow(calcdata80))
calcdata80[ , Z] <- rep(wmeans.80[Z], nrow(calcdata80))
calcdata80 <- as.matrix(calcdata80)
head(calcdata80)

calcdata00 <- (combo.dat[ , c(1, 3:14)])
wmeans.00[X]
wmeans.00[Y]
wmeans.00[Z]
calcdata00[ , X] <- rep(wmeans.00[X], nrow(calcdata00))
calcdata00[ , Y] <- rep(wmeans.00[Y], nrow(calcdata00))
calcdata00[ , Z] <- rep(wmeans.00[Z], nrow(calcdata00))
calcdata00 <- as.matrix(calcdata00)
head(calcdata00)

d_orientation2 <- NULL  

for(i in 1:1999){
  effect_sumissue.80.SP_FL <- exp(calcdata80 %*% bcsimbetas[i, 27:39])/(1 + exp(calcdata80 %*% bcsimbetas[i, 1:13]) + exp(calcdata80%*% bcsimbetas[i, 27:39]) +  exp(calcdata80 %*% bcsimbetas[i, 40:52]))
  effect_sumissue.00.SP_FL <- exp(calcdata00 %*% bcsimbetas[i, 27:39])/(1 + exp(calcdata00 %*% bcsimbetas[i, 1:13]) + exp(calcdata00 %*% bcsimbetas[i, 27:39]) +  exp(calcdata00 %*% bcsimbetas[i, 40:52]))
  diff <- effect_sumissue.00.SP_FL - effect_sumissue.80.SP_FL    
  d_orientation2[i] <- mean(diff)
}

head(d_orientation2)
length(d_orientation2)
mean(d_orientation2)
-0.01616723
sd(d_orientation2)
0.002189141

##p-values:
ate <- mean(d_orientation2)
ate.SE <-  sd(d_orientation2)
tstat <- (ate-0)/ate.SE
pvalue <- 2*(1 - pnorm(abs(tstat), 0, 1))
# 1.523226e-13


##########
##Demographic Profile: (age(4), female(5), black(6), education(3))

X <- 4
Y <- 5
Z <- 6
Q <- 3

calcdata80 <- (combo.dat[ , c(1, 3:14)])
wmeans.80[X]  
wmeans.80[Y]
wmeans.80[Z]
wmeans.80[Q]
calcdata80[ , X] <- rep(wmeans.80[X], nrow(calcdata80))
calcdata80[ , Y] <- rep(wmeans.80[Y], nrow(calcdata80))
calcdata80[ , Z] <- rep(wmeans.80[Z], nrow(calcdata80))
calcdata80[ , Q] <- rep(wmeans.80[Q], nrow(calcdata80))
calcdata80 <- as.matrix(calcdata80)
head(calcdata80)

calcdata00 <- (combo.dat[ , c(1, 3:14)])
wmeans.00[X]
wmeans.00[Y]
wmeans.00[Z]
wmeans.00[Q]
calcdata00[ , X] <- rep(wmeans.00[X], nrow(calcdata00))
calcdata00[ , Y] <- rep(wmeans.00[Y], nrow(calcdata00))
calcdata00[ , Z] <- rep(wmeans.00[Z], nrow(calcdata00))
calcdata00[ , Q] <- rep(wmeans.00[Q], nrow(calcdata00))
calcdata00 <- as.matrix(calcdata00)
head(calcdata00)

d_demo2 <- NULL  

for(i in 1:1999){
  effect_sumissue.80.SP_FL <- exp(calcdata80 %*% bcsimbetas[i, 27:39])/(1 + exp(calcdata80 %*% bcsimbetas[i, 1:13]) + exp(calcdata80%*% bcsimbetas[i, 27:39]) +  exp(calcdata80 %*% bcsimbetas[i, 40:52]))
  effect_sumissue.00.SP_FL <- exp(calcdata00 %*% bcsimbetas[i, 27:39])/(1 + exp(calcdata00 %*% bcsimbetas[i, 1:13]) + exp(calcdata00 %*% bcsimbetas[i, 27:39]) +  exp(calcdata00 %*% bcsimbetas[i, 40:52]))
  diff <- effect_sumissue.00.SP_FL - effect_sumissue.80.SP_FL    
  d_demo2[i] <- mean(diff)
}

mean(d_demo2)
0.001002094
sd(d_demo2)
0.002589598

##p-values:
ate <- mean(d_demo2)
ate.SE <-  sd(d_demo2)
tstat <- (ate-0)/ate.SE
pvalue <- 2*(1 - pnorm(abs(tstat), 0, 1))
# 0.6987791


##########
##Election-Level: (changerdpi(), reelectcamp())

X <- 12
Y <- 13

calcdata80 <- (combo.dat[ , c(1, 3:14)])
wmeans.80[X]  
wmeans.80[Y]
calcdata80[ , X] <- rep(wmeans.80[X], nrow(calcdata80))
calcdata80[ , Y] <- rep(wmeans.80[Y], nrow(calcdata80))
calcdata80 <- as.matrix(calcdata80)
head(calcdata80)

calcdata00 <- (combo.dat[ , c(1, 3:14)])
wmeans.00[X]
wmeans.00[Y]
calcdata00[ , X] <- rep(wmeans.00[X], nrow(calcdata00))
calcdata00[ , Y] <- rep(wmeans.00[Y], nrow(calcdata00))
calcdata00 <- as.matrix(calcdata00)
head(calcdata00)

d_election2 <- NULL  

for(i in 1:1999){
  effect_sumissue.80.SP_FL <- exp(calcdata80 %*% bcsimbetas[i, 27:39])/(1 + exp(calcdata80 %*% bcsimbetas[i, 1:13]) + exp(calcdata80%*% bcsimbetas[i, 27:39]) +  exp(calcdata80 %*% bcsimbetas[i, 40:52]))
  effect_sumissue.00.SP_FL <- exp(calcdata00 %*% bcsimbetas[i, 27:39])/(1 + exp(calcdata00 %*% bcsimbetas[i, 1:13]) + exp(calcdata00 %*% bcsimbetas[i, 27:39]) +  exp(calcdata00 %*% bcsimbetas[i, 40:52]))
  diff <- effect_sumissue.00.SP_FL - effect_sumissue.80.SP_FL    
  d_election2[i] <- mean(diff)
}

mean(d_election2)
-0.002031795
sd(d_election2)
0.002642188

##p-values:
H0: d_election2 = 0
hist(d_election2, breaks = 100)  #sampling distribution of difference in eras.
ate <- mean(d_election2)
ate.SE <-  sd(d_election2)
tstat <- (ate-0)/ate.SE
pvalue <- 2*(1 - pnorm(abs(tstat), 0, 1))
#0.4419041

